Experimental study and simulations of hydrogen cooling effectiveness for aviation PEM fuel cells

Proton exchange membrane fuel cells (PEMFCs) are seen as one possible future means of driving the change towards a zero-emission society. In a civil aircraft, fuel cell systems can have multiple potential benefits, such as reduced noise, lowered emissions and higher fuel economy compared to jet aircraft. For controlling the fuel cell temperature, thermal management systems are required which can be optimized for aircraft applications regarding their weight and reliability. In this work, a simplified and light-weight thermal management system relying on hydrogen cooling is presented and analysed. To investigate the feasibility, a test rig and a three-dimensional, singular channel model in ANSYS Fluent were designed. Fuel cell temperature could be maintained within the set threshold in the model and the test rig, thus showing that controlling the fuel cell temperature via the hydrogen reactant flow is a viable alternative thermal management system. Results from the model indicate that both the hydrogen mass flow and hydrogen inlet temperature should be used to control the fuel cell temperature. Furthermore, operating the fuel cell at medium to low current densities is favourable for hydrogen cooling. Future studies will explore alternate flow field designs to facilitate thermal management system relying on hydrogen.

www.nature.com/scientificreports/heat exchanger onboard fuel cell powered commercial aircraft and found a drag reduction of 23% of the phase change cooling heat exchanger.A recent European patent application by Gao and Kösters 12 presented the idea of hydrogen cooling for fuel cells in aviation, where the hydrogen reactant flow is used both as fuel and coolant, leading to a possible reduction in size and mass of the thermal management system.
Furthermore, system reliability is of great importance in the aviation sector and is often achieved through system redundancy 13 .Currently, a minimum redundancy of 2.25 is required for a fuel cell-powered airliner to pass the certification based on the state-of-the-art fuel cell reliability index in literature 14 .This however increases the total system weight and limits the payload of the aircraft drastically, making them an uncompetitive solution for clean aviation.Hereon, reducing the total number of the system components, i.e., increasing the system simplicity is seen as an effective direction in enhancing the system reliability.Besides, this is accompanied by a weight reduction of the system, leading to an increased power to weight ratio.
In this work, a simplified and thus potentially lightweight thermal management system relying on hydrogen cooling is presented and analyzed regarding its feasibility in a proton-exchange membrane fuel cells (PEMFCs).This is done experimentally and via a numerical model in ANSYS Fluent.

Experimental set-up
All experiments were performed on a PEM fuel cell short stack from ZBT GmbH operated with hydrogen as fuel and oxygen as oxidizer.The stack is water-cooled, consists of six cells and provides electrical power of up to 100 W at 3.6 V.The balance of plant of the fuel cell stack is done by the SMART2 Fuel Cell Test System by WonATech.The SMART2 test rig controls the flow rate of the fuel and oxidizer gasses, records the current and voltage of the fuel cell and measures the pressure and temperature of the fuel and oxidizer gasses.The given accuracy of the gas supply unit is ±1% at a repeatability of ±2% of the full scale.Additionally, the gas streams can be conditioned by heating, humidifying and pressurizing.The fuel cell stack is cooled by the RT 180 thermostat from Lauda through deionized water.The SMART2 test rig is operated through a graphical user interface where the abovementioned parameters can be controlled.Thermal imaging was used to monitor the temperature of the system through a FLIR DM285 multimeter.In the used temperature range, the FLIR DM285 thermal imaging sensor (FLIR Lepton microbolometer) had an accuracy of ±5 K.
The mass flow in the recirculation loop was measured with a SmartTrek 100 mass flow meter from Sierra.The supplier of the mass flow meter states the accuracy as ±1% and the repeatability as ±0.2% of the full scale.The hydrogen is recirculated in the loop by the MH 0018 A hydrogen recirculation pump from Busch Vacuum Solutions.The analogue motor speed control of the hydrogen pump was done by a VSP 2206 from Voltcraft, while the power was supplied by an EA-PS 3040-40 C from EA Elektro Automatik.The hydrogen was cooled through a plate heat exchanger from WilTec.The pipes and fittings to connect the auxiliary equipment were manufactured by Swagelok.A schematic representation of the experimental set-up is shown in Fig. 1.
The hydrogen pump was attached to a wooden board which was supported by foam dampeners for electrical insulation.Furthermore, the piping going to and from the pump was realized with perfluoroalkoxy (PFA) tubes.The required electrical insulation by the supplier of Z Ohm ≥ 50 k was achieved with this set-up.The components of the set-up where hydrogen leakage could occur were fitted inside a ventilated space for safety reasons.A casing of foam insulation was constructed for the fuel cell in order to reduce cooling effects from the ambience and the ventilation system.

Numerical model set-up
In a fuel cell, a wide variety of different physical phenomena occur which are interdependent.In the following, these physical phenomena will be presented with the respective governing equations which are used in the single channel model realized in ANSYS Fluent.

Electrochemistry modeling
The electrochemistry model adopted in ANSYS Fluent is based on work of Kulikovsky et al., Mazumder and Cole, and Um et al. [15][16][17] .The main task of the model is to compute the rate of reactions at the anode and cathode.This is done via the surface overpotential, which is the difference between the phase potential of the solid and the phase potential of the ionomer.The resulting electron transport through the solid conductive materials and the protonic transport through the ionomer is modeled as follows.
Here κ is the conductivity, V the electric potential and j the volumetric transfer current.The subscripts ion and sol denote the values for the ionomer and the solid respectively.∇ is the Nabla operator.The proton conductivity of the ionomer κ ion is modeled according to Springer et al. 18 as follows.
where T is the temperature and l M the water content of the membrane.For the solid phase, j sol = +j cat on the cathode side and j sol = −j an for the anode side.For the ionomer phase, j ion = −j cat on the cathode side and j ion = +j an for the anode side.These terms are so-called source terms and are modeled with the simplified Butler-Volmer equation as follows.
Here i ref is the reference exchange current density, ζ the specific active surface area, γ the concentration dependence, a the charge transfer coefficient, R the universal gas constant and F the Faraday constant.The subscripts cat and an denote the values at the cathode and anode respectively.[A] and [C] represent the molar concentration of the species upon which the reaction rates depend for the anode and cathode respectively.The local surface overpotential η quantifies the voltage loss.On the anode side, it can be modeled over the difference between the solid and ionomer potentials.At the cathode side, the open-circuit voltage V oc is further subtracted from the difference to account for the gain in electrical potential.
For the ionomer phase potential V ion there is a zero flux boundary condition as no ionic current leaves the fuel cell through any external boundary.The solid phase potential V sol also has zero flux boundary conditions for all external boundaries that are not in contact with the outside electrical circuit.For these contact zones, galvanostatic boundary conditions were chosen.

Current and mass conservation
The chemical reaction on the anode and cathode side occurs at the so-called triple-phase-boundary of catalyst, ionomer and species.In the triple-phase boundary at the catalysts on the anode and cathode sides, the creation of the different species is modeled as follows.
Here S is the volumetric source term and M the molar mass of hydrogen ( H 2 ), oxygen ( O 2 ) and water ( H 2 O ) respectively.S H2 and S O2 are negative as they are being consumed at the triple-phase boundary while S H2O is www.nature.com/scientificreports/positive as it is being created.As the total electrical current which is generated at the anode and cathode is equal, the current conversation can be stated as follows.

Heat source
The heat generated in the fuel cell stems from the enthalpy change due to the electrochemical reactions h react , and the voltage losses which are transformed to heat.The voltage losses can be differentiated into activation losses and ohmic losses.Furthermore, the enthalpy change caused by the condensation and vaporization of water h L must be considered.Thus, the heat source term S h can be expressed as follows.
Here I is the electrical current and Z ohm the ohmic resistivity.The heat transfer coefficient α is calculated based on the combined convective and radiative heat flux q , the wall temperature T wall and a predefined refer- ence temperature T ref as follows.

Liquid water formation and transport
The liquid formation and transport model is based on work from 19,20 .In this approach, the conservation equation for the volumetric water saturation s governs the water formation and transport.
Here r H2O is the condensation rate for water, u the velocity, the porosity and ρ the density.The subscript l denotes values for liquid water.The condensation rate is expressed as follows.
Here c r is the condensation rate constant, p H2O the water vapor pressure and p s the saturation pressure.For the porous zones, the convective term in Eq. ( 14) is replaced by a capillary diffusion term.
Here t is the time, K the absolute permeability, µ is the viscosity and p cap the capillary pressure.The capil- lary pressure is computed through the Leverett function and thus is dependent on the water saturation and the wetting phase 21 .
Here θ is the contact angle and Ŵ the surface tension.The effective diffusivity D j eff of the respective gas phase j is modelled via the full multicomponent diffusion method with corrections for the tortuosity of the porous media.
Here D ij is the binary mass diffusion coefficient which is obtained through the Maxwell-Stefan equations for diffusive mass flux which lead to generalized Fick's law diffusion coefficients 22 .
The model was validated through experimental data from Iranzo et al. 23 .As material and geometric parameters of the ZBT fuel cell stack used in the experiments could not be obtained, the given fuel cell parameters from Iranzo et al. were used for the initial validation of the model.Consequently, the operational parameters were also chosen in compliance with the experimental data from Iranzo et al.An overview of the model parameters for the bipolar plates, gas diffusion layers (GDL), catalyst layers and membrane can be found in Table 1.A mesh sensitivity study was conducted and a mesh of 375,000 elements was chosen to minimize the error while still maintaining an adequate model solving time.
The results of the ANSYS fuel cell model compared to the experimental data of Iranzo et al. can be seen in Fig. 2. The percentage wise error reaches a maximum at the end of the analyzed interval of 6.5%.This is satisfactory, especially since the fuel cell system may primarily be operated at lower current densities ranges which are better suited for aircraft applications 2 .Additionally, the polarization curve of the ZBT stack is shown in Fig. 2 which is used in the experimental part of this study.(11)   an j an dV = cat j cat dV

Results and discussion
In the following, the experimental and simulative results will be presented to evaluate the feasibility to achieve sufficient cooling of a PEM fuel cell through the hydrogen reactant stream.First, the overall feasibility will be presented and discussed in the experimental study, followed by an evaluation of beneficial operating parameters through a sensitivity study.

Experimental study
In the following, the results of the parametric, experimental study will be briefly presented.The impact of the hydrogen cooling on the fuel cell stack temperature can be seen in Fig. 3, where experimental and numerical results are combined.The fuel cell was operated at I FC = 5 A , p op = 2 bar , O2 = 2 and H2 = 1.2 prior to initia- tion of hydrogen cooling.The inlet gasses were fully humidified.After stopping the initial water coolant flow, the fuel cell temperature rose from T FC = 42 • C up to T FC = 50 • C over the course of 60 min.It is notable that the fuel cell stack's temperature did not increase consistently over the 60-min interval.While the reaction kinetics inside the fuel cell are temperature dependent, that does not fully explain the occurring plateaus marked in Fig. 3.It is more likely, that an interaction with other systems such as the residual coolant caused an additional but limited cooling effect.After 60 min, the threshold temperature of T FC = 50 • C was reached and hydrogen cooling was initiated.The threshold was chosen so that there is enough buffer to critical membrane temperatures in case hydrogen cooling would not work as anticipated.Initiation of the hydrogen cooling caused a rapid decline of the overall stack temperature and reduces the temperature to T FC = 40 • C over an interval of 30 min.In this experimental set-up, the cooling time of the fuel cell was protracted due to the thermal mass of the residual cooling water in the system.The hydrogen pump was operated at 25% load which resulted in a stoichiometric ratio of H2 = 65.6 .The hydrogen which did not react within the fuel cell was recirculated by said hydrogen pump back to the hydrogen inlet.It can further be seen that the numerical results accurately follow the experimental temperature curve with the maximum temperature difference being below T = 0.5 K.

Numerical model study
In the following, the most representative results of the numerical model study are presented.The entirety of the results of the parametric study can be found in Supplementary Data S1.

Variation of hydrogen stoichiometric ratio
The dependencies of the mean hydrogen channel temperature T cha , the mean hydrogen temperature T H2 and the maximum membrane temperature T M,max as well as the mean heat transfer coefficient of the hydrogen channel α cha to the hydrogen inlet stoichiometry ratio can be seen in Fig. 4a.All temperature values, especially the maximum membrane temperature, are throughout the chosen interval within the operating temperature regions of common PEM fuel cells.Furthermore, the curves of all three temperatures decrease over the interval and flattens out for greater hydrogen stoichiometry ratios as it approaches the inlet hydrogen flow temperature.The mean heat transfer coefficient of the hydrogen channel follows an opposite trend and shows to be sensitive to the hydrogen stoichiometry ratio.The difference between the mean hydrogen temperature and the maximum membrane temperature increases over the interval from T = 2.02 K up to T = 5.66K.
The local distribution of the temperature along the fuel cell channel can be seen for H2 = 80 and H2 = 440 in Fig. 5.The hydrogen flow enters at the upper right corner while the air flow enters at the lower-left corner.It can be seen that the temperature along the fuel cell channel varies strongly with varying hydrogen stoichiometry ratios and is more uniformly distributed for H2 = 80 .This can be quantified by the standard deviation of the mean channel temperature which is minimal for H2 = 80 with σ Tcha = 1.66 K and reaches its maximum at H2 = 200 with σ Tcha = 2.56 K and decreases down to σ Tcha = 2.08 K for H2 = 440.A comparison of the temperature distributions for H2 = 80 and H2 = 440 across the channel at half-length can be seen in Fig. 6.For better comparison, the temperature difference based on the relative temperature of the cathode bipolar plate was plotted.The regions of the fuel cell geometry are marked in this frontal view of the geometry.It can be seen that for H2 = 440 the overall temperature difference between the anode and cathode bipolar plates is more pronounced.Furthermore, the temperature gradient across the channel and the gas diffusion layers is also higher.The maximum temperature in both cases is reached in the cathode channel.
The required blower power is exceptionally sensible to changes of the stoichiometry.The parasitic power consumption of the blower can be displayed as percentage respective to the power produced by the fuel cell channel.For a stoichiometry of H2 = 40 the relative parasitic power amounts to 0.006% but increases to 0.628% for H2 = 440.

Variation of hydrogen inlet temperature
In Fig. 7, a representation of the mean air relative humidity RH air , air relative humidity at outlet RH air,out , mean hydrogen relative humidity RH H2 and hydrogen relative humidity at outlet RH H2,out as a function of the hydrogen inlet temperature T H2,in is presented.It is notable, that the mean hydrogen relative humidity and the outlet hydro- gen relative humidity follow approximately the same upward trend over the chosen interval.Simultaneously, the mean air relative humidity as well as the air outlet relative humidity decrease with increasing values of T H2,in .Especially the mean air relative humidity is relatively high for lower values of the hydrogen inlet temperature with a maximum of RH air = 185% for T H2,in = 273.25 K .It then decreases sharply until being equal to the air outlet relative humidity at around RH air = 101% for T H2,in = 323.15K .This trend can be explained by the cool- ing effect of the hydrogen on the air flow as a decrease in the air temperature causes the relative humidity to rise.red), mean hydrogen temperature T H2 (light blue), maximum membrane temperature T M,max (yellow) and mean heat transfer coefficient of the hydrogen channel α cha (dark green) for varying hydrogen stoichiometry ratios H2 (a) and varying current density i FC (b).www.nature.com/scientificreports/

Variation of fuel cell current density
In Fig. 4b the mean channel temperature T cha , the mean hydrogen temperature T H2 , the maximum membrane temperature T M,max and the mean heat transfer coefficient of the hydrogen channel α cha as a function of the current density i FC can be seen.With increasing current density, the maximum temperature of the membrane increases up to a maximum of T M,max = 316.4K for i FC = 0.6 A cm −2 .This can be explained by the rising heat production due to the increase in generated power and decrease in fuel cell efficiency.Consequently, the mean hydrogen temperature and mean hydrogen channel temperature increase over the set interval and reach at i FC = 0.6 A cm −2 maximum values of T cha = 309.2K and T H2 = 306.9K.The initial decrease of the mean hydrogen temperature might stem from the increased hydrogen mass flow as fuel consumption rises due to increased current density and the hydrogen stoichiometry ratio H2 = 360 was kept constant.The stoichiometric hydrogen factor was chosen based on preliminary simulations, as a higher heat production for greater current densities was anticipated.For greater current densities, the increasing generated heat seems to surpass the increase in cooling capacity due to a higher hydrogen mass flow, leading to an increase in overall temperature.This explanation is supported by the steady increase of the mean heat transfer coefficient Red: i FC = 0.1 A cm −2 at H2 = 360 , light blue: i FC = 0.6 A cm −2 at H2 = 360 , yellow: i FC = 0.4 A cm −2 at H2 = 80 , dark green: i FC = 0.4 A cm −2 at H2 = 440 .For visibility reasons, only a subset of the data is highlighted as markers.α cha over the set interval as the heat transfer is facilitated due to the higher mass stream and the increased tem- perature difference between the hydrogen flow and the membrane.
The temperature distribution across the channel at medium length for current densities of i FC = 0.1 A cm −2 and i FC = 0.6 A cm −2 is shown in Fig. 6.It is visible that for i FC = 0.1 A cm −2 the temperature difference across the channel is with T = 0.98 K below one Kelvin, whereas for i FC = 0.6 A cm −2 the temperature difference totals T = 6.57K .It has to be noted though that this maximum temperature difference is reached when com- paring the hot cathode bipolar plate to the lowest temperature in the hydrogen flow.When it comes to safety concerns however, the maximum temperature difference of the fuel cell components is of greater importance.As can be seen, the maximum temperature difference of the fuel cell components is between the cathode and anode bipolar plates.The temperature difference of the bipolar plates is T = 0.32 K and T = 2.72 K for i FC = 0.1 A cm −2 and i FC = 0.6 A cm −2 respectively.This increase in temperature difference for higher current densities can be attributed to poor thermal conductivity in the MEA and GDL of the fuel cell.

Variation of fuel cell channel geometry
Besides the straight channel geometry presented above, two additional fuel cell channel geometries were studied to quantify the pressure drop along the geometry.A L-and U-shape were chosen to reflect critical geometric aspects of a parallel and serpentine flow field respectively.The three geometries were evaluated for i FC = 0.4 A cm −2 , H2 = 280 p op = 4 bar , T H2,in = 30 • C and fully humidified inlet gasses.The pressure gradi- ent p over the channel for the three studied geometries I-, L-, U-shape for identical operation parameters and channel length was found to be 43.49Pa , 47.69 Pa and 54.68 Pa respectively.The pressure gradient increased from the I-shape to the L-shape by 9.66% and for the I-shape to the U-shape by 25.73%.
The results indicate that a parallel flow field is advantageous for hydrogen cooling as the pressure drop is significantly reduced.Furthermore, the mean channel length of a parallel flow field is smaller than the mean channel length of a serpentine flow field which facilitates the hydrogen cooling as temperature gradients along the respective channel are reduced.Common flow fields are not designed for the high stoichiometry ratios required for hydrogen cooling.Thus, the efficiency of hydrogen cooling can most likely be increased by adapting and optimizing the flow field.

Conclusions and outlook
In order to investigate the feasibility of a solely hydrogen cooled proton-exchange membrane fuel cell, a test rig was constructed and a three-dimensional, singular channel fuel cell model was created in ANSYS Fluent.The model was used to perform a parametric study of the operational parameters to find a potential optimal operation point.
Adequate cooling of the fuel cell was possible while staying within the threshold of the defined safety parameters.The results indicate that varying the stoichiometry ratio of the hydrogen flow and the hydrogen inlet temperature is well suited for controlling the fuel cell temperature.However, especially for longer channel geometries, these two parameters must be handled with great care as both come with respective challenges.On the one hand, high stoichiometry lead to a fast increase of the necessary blower power.On the other hand, changing the inlet temperature of the hydrogen to the lower regions of the interval leads to undesired high-temperature gradients along the channel.This effect might be less pronounced when analyzing a complete flow field instead of a singular straight channel.Furthermore, low inlet hydrogen temperatures require an oversaturated hydrogen flow to prevent dry conditions.This might prove harmful however when analyzing a complete flow field as channels may be flooded.
Although cooling of the fuel cell model via hydrogen cooling was feasible for all chosen current densities, care must be taken when choosing the operating point.As more heat is produced for greater current densities, the cooling of the fuel cell also becomes more challenging.Depending on the application of the fuel cell, lower current densities might be favorable to minimize the overall temperature gradient and the blower power.In addition, it was shown that the thermal conductivity of the fuel cell MEA may be a limiting factor, as the increased temperature gradient across the channel can be largely attributed to poor thermal conductivity of the MEA components.This effect could be drastically reduced by using MEA materials with better thermal conductivity.Moreover, this would also enhance the cooling capabilities of other thermal management systems, as this problem is inherent to the fuel cell system.It is anticipated that a parallel flow field is better suited for hydrogen cooling due to reduced pressure loss and shorter mean channel lengths.
The feasibility of hydrogen cooling was also shown in experiments on the designed test rig for operation at low current densities.After initiating hydrogen cooling, the stack temperature converged rapidly from 50 °C to 40 °C for a hydrogen stoichiometric factor of H2 = 65.6 .This behavior could be replicated in the created numerical model.The results show the viability of hydrogen cooling as an alternative, potentially light-weight thermal management system.Future studies will explore the advantages of alternate flow field design to facilitate hydrogen cooling.Conducting a variation of material parameters regarding the fuel cell stack is also aspired to improve the efficiency of hydrogen cooling for fuel cells.

2 Pt 2 PtFigure 2 .
Figure 2. Validation of the ANSYS fuel cell model (light blue) through comparison with experimental data by Iranzo et al. 23 (yellow).The polarization curve of the ZBT stack is shown in red.

Figure 3 .
Figure 3. Depiction of the fuel cell stack temperature over time in experiment (red) and simulation (blue, dashed).Hydrogen cooling is initiated after 60 min.

Figure 4 .
Figure 4. Depiction of the mean channel temperature T cha (red), mean hydrogen temperature T H2 (light blue), maximum membrane temperature T M,max (yellow) and mean heat transfer coefficient of the hydrogen channel α cha (dark green) for varying hydrogen stoichiometry ratios H2 (a) and varying current density i FC (b).

Figure 5 .
Figure 5. Local distribution of the temperature along the fuel cell channel geometry for H2 = 80 (above) and H2 = 440 (below) at 0.4 A cm −2 .Hydrogen enters the channel at the right upper corner.

Figure 6 .
Figure 6.Comparison of the local temperature difference distribution across the cell at half channel length.Red: i FC = 0.1 A cm −2 at H2 = 360 , light blue: i FC = 0.6 A cm −2 at H2 = 360 , yellow: i FC = 0.4 A cm −2 at H2 = 80 , dark green: i FC = 0.4 A cm −2 at H2 = 440 .For visibility reasons, only a subset of the data is highlighted as markers.

Figure 7 .
Figure 7. Depiction of the mean air relative humidity RH air (red), air relative humidity at outlet RH air,out (light blue), mean hydrogen relative humidity RH H2 (yellow) and hydrogen relative humidity at outlet RH H2,out (dark green) as a function of the hydrogen inlet temperature T H2,in .

Table 1 .
23del parameters used in the ANSYS Fluent model taken from Iranzo et al.23.